plotG(vx,vy,a,b,c)

08783b1b-97e6-4d97-ace5-a0c8eb1afaa1 25 15 12 decimal

&[DATE] &[TIME] - &[FILENAME]

plotG(vx,vy,a,b,c)

vx vy a b c plotG a char ≡ b size ≡ c color ≡ 13 mat n vx length : plot vx 1 el vy 1 el a b c augment : i 2 n range plot plot vx i el vy i el a b c augment stack : for plot 51 line :

data n Decimate collects every 'n' points range 1 data rows n range : k 1 : i range u k el data i 1 el : v k el data i 2 el : k k 1 + : 31 line for u v augment 31 line 31 line :

data

data 0 0.0385 0.2 0.0416 0.4 0.0451 0.6 0.0491 0.8 0.0536 1 0.0588 1.2 0.0648 1.4 0.0716 1.6 0.0796 1.8 0.089 2 0.1 2.2 0.1131 2.4 0.1289 2.6 0.1479 2.8 0.1712 3 0.2 3.2 0.2358 3.4 0.2809 3.6 0.3378 3.8 0.4098 4 0.5 4.2 0.6098 4.4 0.7353 4.6 0.8621 4.8 0.9615 51 5.2 0.9615 5.4 0.8621 5.6 0.7353 5.8 0.6098 6 0.5 6.2 0.4098 6.4 0.3378 6.6 0.2809 6.8 0.2358 7 0.2 7.2 0.1712 7.4 0.1479 7.6 0.1289 7.8 0.1131 8 0.1 8.2 0.089 8.4 0.0796 8.6 0.0716 8.8 0.0648 9 0.0588 9.2 0.0536 9.4 0.0491 9.6 0.0451 9.8 0.0416 10 0.0385 512 mat :

t0 1 time :

data rows 51

X data 1 col : Y data 2 col : 12 mat L X Min : H X Max : 12 mat x c 1 L x ≤ (x H ≤ (& cases : pts X Y . 9 red plotG : 41 sys

pts

===== generalise from here =====

n data rows : N n 1 - : r 2 N range : 31 line

Prepare the cubic splines

Install the constant vector 'br' i r br i el Y i 1 + el Y i el - X i 1 + el X i el - / (Y i el Y i 1 - el - X i el X i 1 - el - / (- eval : br N 1 + el 0 : 21 line for

100000000000000000000000155 mat

Initialise the blank
diagonal matrix M 11 el 1 : M n n el 1 : 21 sys

Part 1. Fill-in the respective l_p_c splines tridiagonal matrices

0. Construct the tridiagonal linear spline matrix [lspline]

i r M i i 1 - el X i el X i 1 - el - 6 / : for

i r M i i el X i 1 + el X i 1 - el - 3 / : for

i r M i i 1 + el X i 1 + el X i el - 6 / : for

1. Construct the tridiagonal parabolic spline matrix [pspline]

M 12 el 1 - :

M N 1 + N el 1 - :

2. Construct the tridiagonal cubic spline matrix [cspline]

M 12 el X 1 el X 3 el - (X 3 el X 2 el - / (:

M N 1 + N el X N el X N 2 - el - (X N 2 - el X N 1 - el - / :

M 13 el X 2 el X 1 el - (X 3 el X 2 el - (/ :

M N 1 + N 1 - el X N 1 - el X N el - (X N 1 - el X N el - (/ :

M M :

the respective spline coefficients s M 1 -^br * eval :

CWD C:\SmathFile CurrentDirectory :

Disable the 'cspline' tridiagonal to collect 'pspline' ... disable 'pspline & cspline' tridiagonal to collect 'lspline.

Spool to file the spline coefficcients s CWD spline wfile 1

s CWD spline rfile :

Part 2. The complimentary cumulative integrate

s 0 0.01173 0.01309 0.01093 0.0182 0.02129 0.01665 0.03212 0.03489 0.03832 0.05182 0.06939 0.07561 0.10816 0.13676 0.16982 0.23398 0.28927 0.37894 0.45995 0.51124 0.43508 0.10343 231 mat

Prepare the pseudo-integrand s X Y fintegrate i 1 N range h i el X i 1 + el X i el - : f i el h i el 2 / Y i el Y i 1 + el + (* h i el (3^24 / s i el s i 1 + el + (* - : 21 line for f 21 line eval :

Σ s X Y fintegrate :

Σ Σ :

Export the cumulative integral Σcumul n X rows 1 - : Sum 1 el 0 : Sum n 1 + el 0 : i 1 n range Sum i 1 + el Sum i el Σ i el + eval : for X Sum augment 51 line :

l_p_c spline cumulative integrate

Σcumul 0 2 1 sys

Part 3. The complimentary scalar dspline

s X Y t dspline t X 2 el ≤ k 1 : t X X rows 1 - el ≥ k X rows 1 - : k 1 : t X k 1 + el > k k 1 + : while 21 line if if c Y k 1 + el Y k el - X k 1 + el X k el - / s k 1 + el 2 s k el * + (X k 1 + el X k el - 6 / (* - : s k 1 + el s k el - X k 1 + el X k el - / (t X k el - (2^2 / * s k el t X k el - (* + c + 31 line eval :

x D s X Y x dspline :

3 D 0.15988

x plotD L x ≤ (x H ≤ (& (x D out spline if :

The scalar dspline

x plotD 0 2 1 sys

2nd derivative of [lpc] spline

X s augment 0 2 1 sys

Part 4. The [lpc] spline interpolation

The while loop form of the spline interpolation function s X Y t interpspline Find the knots that bracket 't' i0 1 : j X rows : j i0 - 1 > i i0 j + 2 / Floor : X i el t > j i : i0 i : if 21 line while Evaluate the spline based on those knots h X j el X i0 el - : a X j el t - h / : b t X i0 el - h / : a Y i0 el * b Y j el * + h 2^6 / a a 2^1 - (* s i0 el * b b 2^1 - (* s j el * + (* + 91 line :

x Spline s X Y x interpspline : x D1 s X Y x dspline : vx X : vy Σcumul 2 col : 12 mat x Σspline s vx vy x interpspline : x D2 s X s x interpspline : x u 1 L x ≤ (x H ≤ (& cases : 61 sys

from dichotomy below var 5.987909 :

x Spline x u * x D1 x u * 0.5 x Σspline * x u * var 2 var 1 - 2 2 mat x D2 x u * 5 1 sys

Sources and credits:

1. Mathsoft for the tridiagonal systems.

2. Robert Adair for the algos:

...InterpSpline(s,X,Y,t)

...dspline(s,X,Y,t)

...fintegrate(s,vx,vy)

3. Smath adapted: Jean Giraud

Robert was "guru" Program Mathcad

]]>

var Spline var D1 var Σspline var Σspline π Σspline - var D2 51 mat 0.5061 0.5066 - 2.1515 1.8554 0.5061 51 mat

2nd derivative of [lpc] spline

X s augment 0 2 1 sys

2nd derivative "scalar'

x D2 x u * 0 2 1 sys

x f x Spline x D2 - :

a 5.75 : b 6.25 : xn b a + (2 / : 31 line

not compatible t Spline t var π int

ε 106 -^:

Dichotomy n 0 : xn f abs ε > xn f a f * 0 ≤ b xn : a xn : if xn b a + (2 / : n n 1 + : 31 line while 31 line

sanity check xn xn f n 31 mat 5.9879 01731 mat

1 time t0 - 860 s *